clear clc;

% 读取心电波形
val = importdata('Ecg.txt');
signal = val(1,1:1800);
fs = 500;
figure;
plot(signal);
title('干净的ECG信号');
xlabel('采样点');
ylabel('幅值db');
grid on;

% 添加高斯白噪声
signal_noise = awgn(signal, 10, 'measured');
figure;
plot(signal_noise);
title('含高斯白噪声的ECG信号');
xlabel('n');
ylabel('幅值db');
grid on;


% %IIR低通滤波器
% %通带截止频率
% wp=0.1*pi;
% %阻带截止频率
% ws=0.16*pi;
% %通带衰减
% Rp=1;
% %阻带衰减
% Rs=15;
% %采样频率
% Fs=500;
% Ts=1/Fs;
% wp1 = 2*tan(wp/2)*Fs;
% ws1 = 2*tan(ws/2)*Fs;
% %选择滤波器最小阶数
% [N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');
% %创建butterworth模拟滤波器
% [Z,P,K]=buttap(N);
% %转为传递函数模型
% [Bap,Aap]=zp2tf(Z,P,K);
% %低通变高通
% [b,a]=lp2lp(Bap,Aap,Wn);
% %双线性法转为数字滤波器
% [bz,az]=bilinear(b,a,Fs);
% %绘制频率响应曲线
% [H,W]=freqz(bz,az);
% figure(3);
% plot(W/pi,20*log10(abs(H)));
% grid;
% xlabel('频率/pi');
% ylabel('dB');
% title('IIR低通滤波器频率响应');


% % FIR低通滤波器
% %通带截止频率
% wp = 0.1*pi;
% % 阻带截止频率
% ws = 0.16*pi;
% %阶数
% wdelta = ws-wp;
% N=ceil(8*pi/wdelta);
% % 截止频率
% wn = (wp+ws)/2;
% %频谱图
% bz=fir1(N-1,wn/pi,hanning(N));
% az = 1;
% [H,W] = freqz(bz,1,512);
% figure;
% plot(W/pi,20*log10(abs(H)));
% grid;
% xlabel('频率/pi');
% ylabel('幅度');
% title('FIR低通滤波器频率响应')



% %添加工频干扰
% av = 100;
% f0 = 50;
% t = 1:length(signal);
% noise = av*cos((2*pi*f0*t/fs));
% signal_noise = noise + signal;
% figure(2);
% plot(signal_noise);
% title('工频干扰的ECG信号');
% xlabel('n');
% ylabel('幅值db');
% grid on;

% % IIR带阻滤波器
% %通带截止频率
% wp1=0.18*pi;
% wp2=0.22*pi;
% %阻带截止频率
% ws1=0.192*pi;
% ws2=0.208*pi;
% %通带衰减
% Rp=1;
% %阻带衰减
% Rs=15;
% %采样频率
% Fs=500;
% Ts=1/Fs;
% wp1 = 2*tan(wp1/2)*Fs;
% wp2 = 2*tan(wp2/2)*Fs;
% ws1 = 2*tan(ws1/2)*Fs;
% ws2 = 2*tan(ws2/2)*Fs;
% % 中心频率和带宽
% w0 = sqrt(wp1*wp2);
% Band = abs(wp1-wp2);
% % 归一化频率
% wp = 1;
% ws = wp*(ws1*Band)/(w0^2-ws1^2);
% %选择滤波器最小阶数
% [N,Wn]=buttord(wp,ws,Rp,Rs,'s');
% %创建butterworth模拟滤波器
% [Z,P,K]=buttap(N);
% %转为传递函数模型
% [Bap,Aap]=zp2tf(Z,P,K);
% %低通变带阻
% [b,a]=lp2bs(Bap,Aap,w0,Band);
% %双线性法转为数字滤波器
% [bz,az]=bilinear(b,a,Fs);
% [H,W]=freqz(bz,az);
% %绘制频率响应曲线
% figure(3);
% plot(W/pi,20*log10(abs(H)));
% grid;
% xlabel('频率/pi');
% ylabel('dB');
% title('IIR带阻滤波器频率响应');


% %添加基线漂移
% n1 = length(signal)/3;
% x1 = zeros(1,n1);
% t = 1:length(signal)-n1;
% x2 = (length(signal)-n1)/2000*(t-1)+1;
% noise = [x1,x2];
% signal_noise = noise + signal;
% figure(2);
% plot(signal_noise);
% title('基线漂移的ECG信号');
% xlabel('n');
% ylabel('幅值db');
% grid on;


% %iir高通滤波器
% %通带截止频率
% wp=0.0028*pi;
% %阻带截止频率
% ws=0.0012*pi;
% %通带衰减
% Rp=1;
% %阻带衰减
% Rs=15;
% %采样频率
% Fs=500;
% Ts=1/Fs;
% wp1 = 2*tan(wp/2)*Fs;
% ws1 = 2*tan(ws/2)*Fs;
% %选择滤波器最小阶数
% [N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');
% %创建butterworth模拟滤波器
% [Z,P,K]=buttap(N);
% %转为传递函数模型
% [Bap,Aap]=zp2tf(Z,P,K);
% %低通变高通
% [b,a]=lp2hp(Bap,Aap,Wn);
% %双线性法转为数字滤波器
% [bz,az]=bilinear(b,a,Fs);
% %绘制频率响应曲线
% [H,W]=freqz(bz,az);
% figure(3);
% plot(W/pi,H);
% grid;
% xlabel('频率/pi');
% ylabel('幅度');
% title('高通滤波器频率响应')




y = filter(bz,az, signal_noise);
figure(4);
plot(y);
title('滤波后的ECG信号');
xlabel('n');
ylabel('幅值db');
grid on;


N = 1800;
signal_noise_XK=fft(signal_noise,N);
%求幅度谱
signal_noise_magXK=abs(signal_noise_XK);
figure(5);
k=0:length(signal_noise_XK)-1;
w = 2*k/N;
stem(w,signal_noise_magXK,'.');
xlabel('频率/pi');
ylabel('Y(K)');
title('滤波前频谱');
grid on;

y_XK=fft(y,N);
%求幅度谱
y_magXK=abs(y_XK);
figure(6);
k=0:length(y_magXK)-1;
w = 2*k/N;
stem(w,y_magXK,'.');
xlabel('频率/pi');
ylabel('Y(K)');
title('滤波后频谱');
grid on;
